00 

o 
o 

Vh 
Oh- 

< 

o 

(N 

o 

B 

a 

o 
o 



> 

m 

^' 

o 

00 

o 



X 



Phase Transitions in Parallel Replication Process 
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The one- dimensional kinetic contact process with parallel update is introduced and studied by 
Monte Carlo simulations. This process is proposed to describe the plant population replication and 
epidemic disease spreading among them. The phase diagram of the model features the line of the 
second order transitions between absorbing and active phases. The numerical results for the critical 
index /3 demonstrate its continuous variation along the transition line accompanied by the variations 
of the structural characteristics of Umiting steady states. We conjecture the non- universality of the 
critical behavior of the model. 
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There has been a steadily growing interest during the 
last two decades or so in the kinetics and phase tran- 
sitions in non-equilibrium systems. The applications of 
those systems range from physics, like, e.g., statistical 
physics, critical phenomena, condensed matter to less 
conventional fields, like biology, ecology or quantitative 
finance, pj, |2|, |3|] The essentially non-equilibrium fea- 
ture of these systems is due to their possibility to irre- 
versibly enter an absorbing (dead, empty) state. The 
mere existence of such a state violates the detailed bal- 
ance. The central problem in studies of non-equilibrium 
systems is the transitions they undergo between various 
active phases and the inactive (absorbing) state. 

The kinetic contact processes which model, e.g., the 
epidemic disease spreading or population replication, 
are the simplest kinetic models exhibiting such non- 
equilibrium phase transitions under variation of their pa- 
rameters. [1, 2., 3.] These variations consist in the change 
of limiting probability distribution for the site occupa- 
tion numbers in replication process, or the infected site 
numbers for epidemic models. Mostly these models pos- 
sess two phases: the absorbing state with the population 
or viruses extinct and the active phase where some sites 
are populated or infected. 

Usually the sequential update formalism is used to 
study these models. It is based on the differential ki- 
netic equation for the probability distribution function. 
In this framework the generic second order transition into 
absorbing state of the directed percolation (DP) univer- 
sality class is established for the majority of such contact 
processes. [1| Another approach is the parallel update 
scheme when the model is represented as a discrete-time 
probabilistic cellular automaton (PC A). [J, |5|, l6| It is 
known that these two approaches can give rather differ- 
ent results [l| and one may reasonably suppose that in 
such case parallel update gives more adequate descrip- 



tion of real kinetics just because in Nature there is no 
queue for the fulfillment of the state transformation pro- 
cesses. So the studies of PCA implementation of replica- 
tion (epidemic) processes may give a more realistic pic- 
ture of them and reveal some new types of critical be- 
havior. 

In this paper we consider a one-dimensional PCA rep- 
resenting the contact process of population replication or 
disease spreading. Imagine a line of plants which main- 
tain population by spreading the seeds to the nearby sites 
so the new plants can grow on them if these sites are 
empty. Let q be the probability of such event for empty 
site having only one neighboring plant. For the case of 
two plants around the empty one we can choose this prob- 
ability r > qiohe r = 1 — (1 — q)'^ = q(2 — q) assuming 
the independence of the two plants' seeds "not spreading 
effects". Some other choices are also possible but here 
we consider only the case r = q{2 — q). Let also p be the 
probability for plant to survive in the one time step. Thus 
we have the PCA on the line of two-state sites (empty 
and full K^ ^i = 0, 1 resp., where i is the intrachain site 
index) with the evolution step probabilities defined for 
the local three-site configurations in Table [D 

TABLE I: Automaton rules 

Si-i,Si,Si+i 111,110,011,010 100,001 101 000 
Probability (5'^ = l) p q q{2 - q) 

Obviously, this model describes also the epidemic pro- 
cess among the plants, since instead of the plant seeds we 
can consider the virus spreading. Note also that this pro- 
cess differs essentially from the one of the cell replication 
considered in Ref.[3]. 

We have performed Monte Carlo simulations of this 
PCA starting from several random initial states for the 



chains of N = 10000 sites for up to 5000 time steps. The 
periodic boundary conditions were implemented. The 
resulting p — q diagram (Fig. [1]) consists of a single tran- 
sition line between absorbing (all occupation numbers 
5'i = in the infinite-time limit) and active (some Si ^ 0) 
phases. The curve has the following approximate analyti- 
cal representation q = 0.66 — 0.913(p — 0.15)^, determined 
from a direct fit. 
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FIG. 1: (Color online) Phase diagram. 

We assume the second order transition between the 
phases, since the concentration of active sites defined as 



m 



1 ^ 



(1) 



appears to be continuous. The examples of such behav- 
ior are presented in Fig. [21 where the (continuous) q- 
dependence of ni is shown for p = 0. 2, 0.5, 0.8. We found 
that the points obtained in simulations follow closely the 
power law rii oc {q — qcY with the indices [3 given in 
Table [m 



TABLE II: Transition points qc and critical indices f3 for dif- 
ferent p values. 



p 


0.2 


0.5 


0.8 


qc 


0.645 


0.535 


0.275 


P 


0.332 


0.286 


0.254 



One can notice the unusual variation the index (3 with q 
instead of its more conventional universality. The index 
varies closely to the directed percolation value jSdp = 
0.276. [l|,l2] A naive explanation of these deviations could 
have been as an artefact of low data precision, while the 
PCA studied belongs actually to the DP class. Indeed, 
there are small variations of ni for different trials (and 
initial states) as well as its fiuctuations in the nominally 
steady state at large times. Not too close to the transition 
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FIG. 2: (Color online) g-dependence of ni for p — 0.2 (x), 0.5 
(o), 0.8 (n). Solid lines show power law dependencies with 
the indices given in Table HIl 



point i\q — qc\ ^ 0.05) these variations of n\ are of the 
order of several percents and could not possibly account 
for the differences in /3 of the order of ten percents. 

Therefore we are lead to the assumption that this repli- 
cation process does not belong to the DP universality 
class, and it is characterized by the variable critical in- 
dices. To date the signs of the non-universal behavior 
are found only in the \d pair contact process with dif- 
fusion [7] but the type of its critical anomalies is still 
debated; see Ref. [8] and references therein. The non- 
universality of some equilibrium spin models (i.e., the 
coupling-dependent critical exponents) is well known. 
For example, it is established for the classical XY model 
[9|], the eight- vertex model solved by Baxter [10], or 
for several two-dimensional Ising models with competing 
nearest-neighbor (nn) and next-nearest-neighbor (nnn) 
interactions, see Refs. [Il|,ll2|,ll3|,ll^ and more references 
therein. More examples and discussions on the equiva- 
lence between the two-dimensional classical spin models 
(Ising, Potts, Ashkin- Teller, XY) and (1+1) quantum 
models (quantum spin chains, Luttinger, Gaussian), their 
critical properties and (non) universality can be found in 
Refs. [i5|,[3. 

From the renormalization group (RG) point of view the 
non-universality of the PCA implies that the RG fiow 
for the equivalent equilibrium spin model [5|, l^] would 
have a manifold of fixed points instead of a unique fixed 
point. A paradigmatic example of such behavior is the 
Kosterlitz-Thouless picture of RG fiow for the XY model, 
LSI which appears also in many (1+1) quantum models. 
|l7j According to Kadanoff and Wegner, non-universality 
can be traced back to the presence of an RG marginal 
operator. |ll5 

To corroborate our conjecture of non-universality we 
point out that the structure of the active state in this 
model undergoes considerable qualitative changes under 



p and q variations at constant rii. The space-time pat- 
terns of the process which has reached a steady state with 
ni ~ 0.5 are presented in Fig. [3]for different p and q. The 
evolution of 100 sites chosen out of 10000 is shown there 
during 100 time steps in the steady state. 
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FIG. 3: (Color online) (a) - p =0.2, q = 0.89, (b) - p =0.5, 
q = 0.66, (c) - p =0.8, ^ = 0.32. Blue sites are empty. The 
time direction is horizontal. 

The pattern difference, quite distinct visually, can be 
assessed numerically. It appears that along with ni, the 
concentration of clusters of adjacent active and inactive 



1 ^ 



N 



(2) 



i=l 



stays almost constant with negligible fluctuations, see 
FiglH In the above equation and throughout, we use 
5(m^n) as a notation for the conventional Kronecker 
delta. 
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FIG. 4: (Color online) Time dependence of ni (o red), ric (□ 
blue), mi (x green), and niii,i (o magenta); (a) p — 0.2, q 
= 0.89, (b) p = 0.8, q = 0.32. 



This parameter, as well as ni, is almost the same in the 
different trials. The values of Uc for three steady states 
in Fig. [3] are presented in Table lllli The concentration 
of active adjacent pairs (11) in steady configurations nn 
is also shown there with 



1 ^ 

11 = T7X]^(^i + 5^+1,2) 



i=l 



Note that an exact relation 



^11 =ni- nc/2 



(3) 



(4) 



TABLE III: Parameters of three steady state configurations 
shown in Fig. [S] 



(p^q) 


(0.2, 0.89) 


(0.5, 0.66) 


(0.8, 0.32) 


ni 


0.499 


0.506 


0.501 


Tie 


0.640 


0.468 


0.354 


nil 


0.179 


0.272 


0.324 


niii,i 


0.016 


0.072 


0.179 



holds. We have also calculated the average values of the 
supposed marginal operator 



^111,1 



1 ^ 



Si-\-i,tSi-i^tSi^t-\-i 



(5) 



This four-spin operator is present in the Hamiltonian of 
the corresponding equilibrium spin model and it can be 
responsible for the non-universality in (1 + 1) dimensions. 
The average nm^i is also nearly constant in steady states 
(see Fig. |4]), trial independent and varies with p and q^ 
as shown in Table IIIIl Its value defines the structural 
characteristics of emergent active state and, probably, 
their indices. Our data also show that nu (and Uc) is 
proportional to rii near the transition points. 

Summary & Discussion: The structural characteristics 
of steady states with nearly the same ni are considerably 
different, which is not only seen quite clearly from Fig. [31 
but it manifests quantitatively via variations of numbers 
of clusters per site ric ^ and the four-spin operator nm^i 
([5]). Along with the proportionality of nu (and Uc) to ni 
this suggests a multi-component order parameter. The 
latter violates one of the requirements of the Janssen- 
Grassberger hypothesis for a model to belong to the DP 
universality class [8]. We conjecture the non- universality 
in this model, most likely due to the marginal four-spin 
term of the Hamiltonian [cf. Eq. ([5])], which manifests 
itself in variations of the critical index /3. 

However to confirm that this PGA does exhibit the 
non-universal critical behavior, the numerical values of 
other critical indices are needed. Although it is a rather 
difficult task for the present model to get them due to 
strong fluctuations which other quantities of interest such 
as correlation lengths exhibit near the transition point. 
Strong fluctuations also hinder determination of the in- 
dex a = 2 — u\\ — u± [1] from the energy-like behavior of 

^111,1 ^ (<7-^c)^"". 

From the analytical side, the non-universality can 
be probed by analysis of the marginal perturbation 
around an integrable fixed point. For instance, the non- 
universality of the eight- vertex model or the Ising with 
competing nn and nnn couplings, can be demonstrated 
by the RG analysis of the flow generated by the nn or 



four-spin interactions which couple two independent (in- 
tegrable) Ising models. Il2|, ll3|, ll5|, ll6| In the present case 
the RG analysis of non- universality, i.e., of the marginal 
perturbations, is more difficult problem, since even the 
(1+1) DP fixed point presumably controlling the univer- 
sality, is not integrable. [8] Note that the e-expansions 
around the DP upper critical dimension dc = 4 are not 
reliable for our case of (i = 1. [18] Glearly, an RG study of 
the field theory corresponding to our model is warranted. 
It can shed more light on the the critical properties of 
the model. In particular, the field theory formulation of 
the present PGA model could help to answer a natural 
question of how this PGA is distinct from others known 
to belong to the DP universality class. pL,] We plan to 
address all these issues in the forthcoming paper. 
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